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Abstract 

We address the common problem of calculating intervals in the presence of 
systematic uncertainties. We aim to investigate several approaches, but here 
describe just a Bayesian technique for setting upper limits. The particular ex- 
ample we study is that of inferring the rate of a Poisson process when there are 
uncertainties on the acceptance and the background. Limit calculating software 
associated with this work is available in the form of C functions. 

1 The problem 

A very common statistical procedure is obtaining a confidence interval for a physics pa- 
rameter of interest, when there are uncertainties in quantities such as the acceptance of 
the detector and/or the analysis procedure, the beam intensity, and the estimated back- 
ground. These are known in statistics as nuisance parameters, or in Particle Physics 
as sources of systematic uncertainty. We assume that estimates of these quantities are 
available from subsidiary measurements.^ A variant of this procedure which is particu- 

"'^ There are other possibilities. Thus it may be that all that is known is that a nuisance parameter 
is contained within a certain range: /xi < /x < yitu; that is not enough information for a Bayesian 
approach. Alternatively the data relevant for the physics and nuisance parameters could be bound up 
in the main measurement, and not require a subsidiary one. 
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larly relevant for Particle Physics is the extraction of an upper limit on the rate of some 
hypothesized process or on a physical parameter, again with systematic uncertainties. 

To specify the problem in more detail, we assume that we are performing a counting 
experiment in which we observe n counts, and that the acceptance has been estimated 
as eo ± and the background as 60 =t CTfc- For a signal rate s, n is Poisson distributed 
with mean se + h. Here e contains factors like the intensity of the accelerator beam(s), 
the running time, and various efficiencies. It is constrained to be non-negative, but can 
be larger than unity. 

We aim to study and compare different approaches for determining confidence in- 
tervals for this problem. In general we are interested in pathologies in these areas: 

• Coverage. This is a measure of how often the limits that we deduce would 
in fact include the true value of the parameter. This requires consideration of 
an ensemble of experiments like the one we actually performed, and hence is 
an essentially frequentist concept. Nevertheless, it can be applied to a Bayesian 
technique. 

Coverage is a property of the technique, and not of the particular limit deduced 
from a given measurement. It can, however, be a function of the true value of 
the parameter, which is in general unknown in a real measurement. 

Undercoverage (i.e. the probability of containing the true value is less than the 
stated confidence level) is regarded by frequentists as a serious defect. Usually 
coverage is required for all possible values of the physical parameter.^ In contrast, 
overcoverage is permissible, but the larger intervals result in less stringent tests 
of models that predict the value of the parameter. For measurements involving 
quantised data (e.g. Poisson counting), most methods have coverage which varies 
with the true value of the parameter of interest, and hence if undercoverage is to 
be avoided, overcoverage is inevitable. 

Frequentist methods by construction will not undercover for any values of the 
parameters. This is not guaranteed for other approaches. For example, even 
though the Bayesian intervals shown here do not undercover, in other problems 
Bayesian 95% credible intervals could even have zero coverage for some values of 
the parameter of interest. [2] It should also be remarked that, although coverage 
is a very important property for frequentists, on its own exact coverage does not 
guarantee that intervals have desirable properties (for many examples, see Refs. 
and 01). 

^The argument is that the parameter is unknown, and so we wish to have coverage, whatever its 
value. This ensures that, if we repeat our specific experiment many times, we should include the 
true value within our confidence ranges in (at least) the stated fraction of cases. This argument may, 
however, be over-cautious. The location of the dips in a coverage plot like that of Fig.noccur at values 
which are not fixed in s, but which depend on the details of our experiment (such as the values of e and 
h). These details vary from experiment to experiment. Thus we could achieve 'no undercoverage for 
the ensemble of experiments measuring the parameter s\ even if the individual coverage plots did fall 
below the nominal coverage occasionally. Thus in some sense 'average coverage' would be sufficient 
(see for example reference PP), although it is hard to quantify the exact meaning of 'average'. It 
should be stated that this is not the accepted position of most High Energy Physics frequentists. 
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• Interval length. This is sometimes used as a criterion of accuracy of intervals, 
in the sense that shorter intervals have less probability of covering false values 
of the parameter of interest. However, one should keep in mind that short in- 
tervals are only desirable if they contain the true value of the parameter. Thus 
empty intervals, which do occur in some frequentist constructions, are generally 
undesirable, even when their construction formally enjoys frequentist coverage. 

Intervals that fill the entire physically allowed range of the parameter of interest 
may also occur in some situations. Examples of this behavior are given in |3I 
and jS]. An experimenter who requests a 68% confidence interval, but receives 
what appears to be a 100% confidence interval instead, may not be satisfied 
with the explanation that he is performing a useful service in helping to keep 
the coverage probability — averaged over his measurement and his competitor's 
measurements — from dropping below 68%. 

• Bayesian credibility. In some situations it may be relevant to calculate the 
Bayesian credibility of an interval, even when the latter was constructed by fre- 
quentist methods. This would of course require one to choose a prior for all 
the unknown parameters. The question is one of plausibility: given the type of 
measurement we are making, the resolution of the apparatus, etc., how likely is 
it that the true value of the parameter we are interested in lies in the calculated 
interval? Does this differ dramatically from the nominal coverage probability of 
the interval? In fact for different values of the observable(s), frequentist ranges 
are very likely to have different credibilities. Some examples of this behavior are 
noted in Ref. 0. 

When calculating the Bayesian credibility of frequentist intervals, "uninforma- 
tive" priors appear advisable. Note that severe interval length pathologies will 
automatically produce a large inconsistency between the nominal coverage and 
the Bayesian credibility of an interval. Except in a handful of very special cases, 
it is not possible to construct an interval scheme that has simultaneously con- 
stant Bayesian credibility and constant frequentist coverage, even if one has total 
freedom in choosing the prior (s). Although it is not at all clear exactly how large 
a level of disagreement is pathological, nevertheless it may be instructive to know 
how severely an interval scheme deviates from constant Bayesian credibility (and 
how sensitive this is to the choice of prior). 

• Bias. In the context of interval selection, this means having a larger coverage 
B{s', Strue) for an incorrect value s' of the parameter than for the true value Strue- 
This requires plots of coverage versus s' for different values of Stme- For upper 
limits, B{s[,StTuc) > B{s2,Strnc) if s[ is less than s'a, so methods are necessarily 
biassed for low s'. Bias thus is not very interesting for upper limits. It will be 
discussed in later notes dealing with two-sided intervals. 

• Transformation under reparametrisation. Intervals that are not 
transformation-respecting can be problematic. For example, it is possible for 
the predicted value of the lifetime of a particle to be contained within the 90% 
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interval determined from the data, but for the corresponding predicted value of 
the decay rate (equal to the reciprocal of the predicted lifetime) to be outside the 
90% interval when the data is analysed by the same procedure, but in terms of 
decay rate. This would result in unwanted ambiguities about the compatibility 
of the data with the prediction. 

• Unphysical ranges. The question here is whether the interval construction 
procedure can be made to respect the physical boundaries of the problem. For 
example, branching fractions should be in the range zero to one, masses should not 
be negative, etc. Statements about the true value of a parameter should respect 
any physical bounds. In contrast, some methods give estimates of parameters 
which can themselves be unphysical, or which include unphysical values when 
the errors are taken into account. We do not recommend truncating ranges of 
estimates of parameters to obey such bounds. Thus the fact that a branching 
fraction is estimated as 1.1±0.2 conveys more information about the experimental 
result than does the statement that it lies in the range 0.9 to 1. 

• Behavior with respect to nuisance parameter. We would normally expect 
that the limits on a physical parameter would tighten as the uncertainty on a 
nuisance parameter decreases; and that as this uncertainty tends to zero, the 
limits should agree with those obtained under the assumption that the "nuisance 
parameter" was exactly known. (Otherwise we could sometimes obtain a tighter 
limit simply by pretending that we knew less about the nuisance parameter than 
in fact is the case.) These desiderata are not always satisfied by non-Bayesian 
methods (see [Zj and [H]). 

Although we are ultimately interested in comparing different approaches to this 
problem, in this note we investigate a Bayesian technique for determining upper limits. 
Our purpose is to spell out in some detail how this approach is used, and to discuss 
some of the properties of the resulting limits in this specific example. We believe that, 
for variants of this problem (e.g. different choice of prior for s; alternative assumptions 
about the information on the nuisance parameters; etc.), the reader could readily adapt 
the techniques described here (and the associated software) to their particular situation. 

We will report on two-sided intervals and also compare with other methods (e.g. 
Cousins-Highland, pure frequentist, profiled frequentist) in later notes. 

2 Reminder of Bayesian approach 

Before dealing with the problem of extracting and studying the limits on s as deduced 
from observing n events from a Poisson distribution with mean se + 6 in the presence 
of an uncertainty on e, we recall the way the Bayesian approach works for the simpler 
problem of a counting experiment with no background and with e exactly known. Then 
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n is Poisson distributed with mean se, and Bayes' Theorem^ 



P{B\C) = P{C\B)P{B)/P{C) 



(1) 



gives 



p{s\n) 



P{n\s)7i{s) 



(2) 



/ P{n\s)7!'{s) ds 



where 7i{s) is the prior probabihty density for s; p{s\n) is the posterior probabihty 
density function (p.d.f.) for s, given the observed n; and P{n\s) is the probabihty of 
observing n, given s. 

We assume a constant prior for s,^ and that P{n\s) is given by the Poisson 



The hmit is now obtained by integrating this posterior p.d.f. for s until we achieve the 
required fraction (3 of the total integral from zero to infinity. If (3 is 90%, the upper 
limit Su is given by 



/? is termed the credible or Bayesian confidence level for the limit. 

For different observed n, the upper limits are shown in the last two columns of 
Table d for 6 = and for 6 = 3 respectively. The Gaussian approximation for the case 
b = 0, n = 20, would yield Su — 20 + 1.28v^20 ~ 25.7, which is roughly comparable to 
the corresponding = 27.0451 of the Table. For 6 = 0, it coincidentally turns out that, 
for this particular example, the Bayesian upper limits are identical with those obtained 
in a frequentist calculation with the Neyman construction and a simple ordering rule 
(see later note on the frequentist approach to this problem). In general this is not so. 
Other priors sometimes used for s are l/i/s [HI or 1/s jTU]. Having a prior peaked at 
smaller values of s in general results in tighter limits for a given observed n. 

If the whole procedure is now repeated with a background b and a flat prior, the 
upper limits not surprisingly decrease for increasing b at fixed n (except for the case 
n = where the limits can trivially be seen to be independent of b). This is not 
inconsistent with the fact that the mean limit for a series of measurements increases 
with b, i.e. experiments with larger expected backgrounds have poorer sensitivity. 

■^We follow the common convention whereby lower case tt's denote prior p.d.f. 's, lower case p's 
denote other p.d.f. 's, upper case II's denote prior probabilities, and upper case P's denote other 
probabilities. Equation is true for probabilities, p.d.f. 's, or mixtures depending on whether B 
and/or C are discrete or continuous variables. 

■^This is an assumption, not a necessity, and is in some ways unsatisfactory. (It is implausible, 
cannot be normalised, and creates divergences for the posterior if used with a (truncated) Gaussian 
prior for the acceptance e.) 

^It turns out that the sum over n of the discrete distribution (j^J and the integral over s of the 
continuous distribution are both equal to unity. This means that the probability P{n\s) and the 
probability density p{s\n) are correctly normalised. 
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1 ±0 


n 


b = 


1 


2 


3 


4 


5 


6 


7 


8 




3 





2.3531 


2.3531 


2.3531 


2.3531 


2.3531 


2.3531 


2.3531 


2.3531 


2.3531 




2.3026 


1 


3.9868 


3.3470 


3.0620 


2.9019 


2.8000 


2.7297 


2.6783 


2.6391 


2.6083 


3.8897 


2.8389 


2 


5.4669 


4.5520 


3.9676 


3.6026 


3.3623 


3.1953 


3.0736 


2.9816 


2.9099 


5.3223 


3.5228 


3 


6.8745 


5.8618 


5.0463 


4.4644 


4.0571 


3.7666 


3.5534 


3.3922 


3.2671 


6.6808 


4.3624 


4 


8.2380 


7.1964 


6.2451 


5.4751 


4.8914 


4.4569 


4.1313 


3.8832 


3.6904 


7.9936 


5.3447 


5 


9.5714 


8.5213 


7.5063 


6.6022 


5.8579 


5.2719 


4.8180 


4.4660 


4.1904 


9.2747 


6.4371 


6 


10.8826 


9.8288 


8.7885 


7.8047 


6.9344 


6.2066 


5.6184 


5.1499 


4.7772 


10.5321 


7.5993 


7 


12.1766 


11.1203 


10.0703 


9.0460 


8.0904 


7.2450 


6.5289 


5.9387 


5.4586 


11.7709 


8.7958 


8 


13.4570 


12.3984 


11.3441 


10.3014 


9.2952 


8.3635 


7.5374 


6.8300 


6.2380 


12.9947 


10.0030 


9 


14.7261 


13.6655 


12.6085 


11.5575 


10.5247 


9.5365 


8.6250 


7.8142 


7.1136 


14.2060 


11.2085 


10 


15.9858 


14.9233 


13.8641 


12.8090 


11.7630 


10.7415 


9.7701 


8.8758 


8.0775 


15.4066 


12.4073 


11 


17.2375 


16.1732 


15.1121 


14.0542 


13.0017 


11.9621 


10.9525 


9.9966 


9.1170 


16.5981 


13.5983 


12 


18.4823 


17.4163 


16.3533 


15.2934 


14.2371 


13.1881 


12.1560 


11.1582 


10.2162 


17.7816 


14.7816 


13 


19.7210 


18.6535 


17.5887 


16.5269 


15.4682 


14.4139 


13.3692 


12.3452 


11.3588 


18.9580 


15.9580 


14 


20.9545 


19.8854 


18.8191 


17.7554 


16.6946 


15.6373 


14.5856 


13.5459 


12.5302 


20.1280 


17.1280 


15 


22.1832 


21.1127 


20.0448 


18.9795 


17.9169 


16.8572 


15.8014 


14.7528 


13.7187 


21.2924 


18.2924 


16 


23.4078 


22.3359 


21.2665 


20.1996 


19.1353 


18.0737 


17.0151 


15.9612 


14.9161 


22.4516 


19.4516 


17 


24.6286 


23.5553 


22.4845 


21.4161 


20.3502 


19.2868 


18.2261 


17.1689 


16.1172 


23.6061 


20.6061 


18 


25.8459 


24.7714 


23.6992 


22.6294 


21.5619 


20.4969 


19.4344 


18.3747 


17.3189 


24.7563 


21.7563 


19 


27.0601 


25.9844 


24.9109 


23.8397 


22.7708 


21.7042 


20.6400 


19.5784 


18.5198 


25.9025 


22.9025 


20 


28.2715 


27.1946 


26.1199 


25.0474 


23.9770 


22.9090 


21.8432 


20.7799 


19.7191 


27.0451 


24.0451 



Table 1: Upper 90% limits for n observed events with b background and e = 1.0 ± 0.1 
[k = 100 and m = 99, as defined in section IT^ . Also shown are limits for 6 = and 
6 = 3 with fixed e = 1. 

2.1 Coverage 

Next we can investigate the frequentist coverage C(strue)^ of this Bayesian approach. 
That is, we can ask what the probability is, for a given value of Strue, of our upper limit 
being larger than Strue, and hence being consistent with it. This is equivalent to adding 
up the Poisson probabilities of eqn. (jS)) for those values of n for which Su_{n) > Stme i-e. 

C{Strne) = e-'^^' {StrneeT / ul (6) 

relevant n 

As Strue increases through any of the values of of the last two columns of Table [U 
the coverage drops sharply. For example, for the case of zero background and efficiency 
known to be unity, the 90% Bayesian upper limits will include Strue = 3.8896 for n = 1 
or larger. But Strue = 3.8898 is no longer below the upper limit for n = 1. Thus 
one term drops out of the summation of eqn. (jU)) for the calculation of the coverage at 
■Strue = 3.8898, while the remaining terms change but little for the small change in Stme] 
this produces the abrupt fall in coverage. The coverage is plotted in Fig. [TJ where the 
drop at Strue = 3.8897 can be seen. 

The calculation of C(struc) can be done as follows: The identity 

-e-- for /(.) = e-Eir 

^This is the coverage at s = Struo when the Poisson variable is generated with s 
differs from B{s' , struo) where the coverage is checked at s = s' when the generation value is stn: 



Ex" X - a;*- 
,A;=0 fc=0 
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allows us to write (integrating ~f'{x)) 

r^mn) ds = = 1 - e--'--^ ^ ^-^^^ (8) 

Jo Jo n. k. 

From this, it follows that "relevant n" is equivalent to any one of these inequalities: 

Sn{n) > Struo ^ p(s|n) ds > p(s|n) ^ /3 > 1 - e"^'™^ (ftrue£)^ 

and our expression for the coverage becomes 

C(st™e)-1-E'^""'""^^^^ (10) 

•■^ — ' n! 

ri=0 

where means "sum until the next term would cause the sum to exceed 1 — /3" . This 
result proves that C(struo) ^ /? for all values of struo in this simple example. 

It is seen that the coverage starts at 100% for small Strue- This is because even for 
n = the Bayesian upper limit will include Stme, and this is even more so for larger n. 

Bayesian methods can be shown to achieve average coverage. By this we mean that 
when the coverage is averaged over the parameter s, weighted by the prior in s, the 
result will agree with the nominal value (3, i.e. 

fCis) nis)ds (11) 
J 7i{s) as 

A proof of this theorem is given in the second appendix, section [7| of this note. 

For a constant prior, the region at large s tends to dominate the average, while in 
general we will be interested in the coverage at small s. Thus the "average coverage" 
result is of academic rather than practical interest, especially for the case of a flat prior. 
Indeed it is possible to have a situation where the average coverage is, say, 90%, while 
the coverage as a function of s is always larger than or equal to 90%. 



3 The actual problem 

Our actual problem differs from the simple case of Section |21 in that 

(a) we have a background b, assumed for the time being to be accurately known; and 

(b) we have an acceptance e estimated in a subsidiary experiment as eo ± cr^. 

What we are going to do is to use a multidimensional version of Bayes' Theorem to 
express p{s, e\n) in terms of P{n\s, e) and the priors for s and e. The relationship is''' 

p^s, e\n) = Pin\s,e)7r{s)n{e) ^^^^ 

JJ P{n\s,e)TT{s)7i{e) ds de 



^For the case where the probabilities have a frequency ratio interpretation, this is seen from the 
mathematical identities 

P{X and Y and Z) - ""j^(^)""'' "Ir^f^^ = P{X,Y\Z) P{Z) and 
P{X and Y and Z) - ""^^if^^^Y) ''^''nTJ = -^(^1^'^) ^(^'^)- 

So with X, Y and Z identified with s, e and n respectively, and with the prior for s and e 
factorising into two separate priors for s and for e, we obtain p[s, e\n) P{n) = P(n\s, e) 7r(s) 7r(e). 
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To obtain the posterior p.d.f. for s, we now integrate this over e: 



j'OO 

p{s\n) = I p(s,e|n)(ie, (13) 
Jq 

and finally we use this to set a limit on s as in eqn. (jS)). 

The coverage for this procedure needs to be calculated as a function of Strue and 
etrue- The average coverage theorem of the previous section must be generalized to 

//C(g,e)7r(g)7r(e) dsde _ 

Jj7r{s)7i{e)dsde ^ ^ ' 



3.1 Priors 

To implement the above procedure we need priors for s and e. As in the simple 
example of Section |21 for simplicity we assume that the prior for s is constant. It will 
be interesting to look at the way the properties of this method change as other priors 
for s are used. 

We assume that the prior for e is extracted from some subsidiary measurement 
eo ± cr^. We do not assume that this implies that our belief about etrue is represented 
by a Gaussian distribution centred on eo, as this would give trouble with the lower 
end of the Gaussian extending to negative e. Instead, we specify some particular form 
of the subsidiary experiment that provides information about e, and then assume that 
a Bayesian analysis of this yields a posterior p.d.f. for e. Slightly confusingly, this 
posterior from the subsidiary experiment is used as the prior for the application of 
Bayes' Theorem to extract the limit on s (see eqns. (|12p and (|13p). 



3.2 The subsidiary measurement 

Somewhat arbitrarily, we assume that, for a true acceptance etrue, the probability for 
the measured value eo in the subsidiary experiment is given by a Poisson distribution 

P(eoletrue) = e-'^^'-^'^eTUml (15) 

where eo = (m + l)//t, erf = (m + and k is a scaling constant^. We interpret 

this as the probability for eo. This is discrete because the observable m is discrete, but 
the allowed values become closely spaced for large k. For small a^/e (i.e. for large m), 
these probabilities approximate to a narrow Gaussian (see Fig. 12)). 

Given our choice of probability in eqn. (|15p. the likelihood for the parameter e, 
given measured eo, is 

C{e\eo) = e~^'K"'e"'/m\ (16) 

This is the same function of e and eo as eqn. ()15|) . but now m is regarded as fixed, and 
e is the variable. The likelihood is a continuous function of e. It is compared with a 
Gaussian in Fig. |21 

^Here we define eo and as the mean and variance of the posterior p.d.f. of eqn. (|17|l . 
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Finally in the Bayes approach, with the choice of a constant prior for e, the posterior 
probability density for e after our subsidiary measurement is 

p{e\m) oc e-"^/€'"e™/m! (17) 

which is obtained by multiplying the right-hand side of eqn. ()16|1 by unity. This pos- 
terior probability density for e will be used as our prior for e in the next step of 
deducing the limit for s. 



4 Results 

The details of the necessary analytical calculations^ are presented in the Appendix of 
this note. In this section we investigate the behavior of the Bayesian limits in this 
example, especially the shape of the frequentist coverage probability as a function of 



4.1 Shape of the posterior 

The posterior p.d.f. for s has the form 



p{s\b, n)ds oc 



de 



Ids 



where the likelihood, the prior for e, the (constant) prior for s, and the marginalization 
integral over e are all prominently displayed. 

The posterior probability density for s gives the complete summary of the outcome 
of the measurement in the Bayesian approach. It is therefore important to understand 
its shape before proceeding to use it to compute a limit (or extract a central value and 
error-bars) . 

Figure m illustrates the shape of the posterior for s (i.e. marginalized over e) in the 
case of a nominal 10% uncertainty on e, and an expectation of 3 background events. 
Plots are shown for 1, 3, 5, and 10 observed events. The posterior evolves gracefully 
from being strongly peaked at s = to a roughly Gaussian shape that excludes the 
neighborhood near s = with high probability. Technically, the posterior would be 
described as a mixture of n + 1 Beta distributions of the 2nd kind^°, giving it a tail at 
high s that is heavier than that of a Gaussian. 



4.2 Upper Umits 

In this note, our main goal is to obtain a Bayesian upper limit from our observation 
of n events. It is by integrating the posterior p.d.f. out to s = Su that an upper limit is 
calculated: a /? = 90% upper limit is defined so that the integral of the posterior from 

^This example can be handled analytically. More complicated cases might require numerical inte- 
gration, which can be done via numerical quadrature or Monte Carlo methods. 

^°The 2nd Beta distribution is also known as "Beta'" (i.e. "Beta prime" ) , "inverted Beta" , "Pearson 
Type VI" , "Variance-Ratio" , "Gamma-Gamma" , "F" , "Snedecor" , "Fisher-Snedecor" .... 
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s = to s = Su is 0.9. The probability (in the Bayesian sense) of Strue < Su is then 
exactly (3. 

TableQshows the upper limits {f3 = 0.9) for n = 0-20 observed events with b = 0-8 
and e = 1.0 ± 0.1. (Integer values of b are chosen for illustration purposes only; b can, 
of course, take any real value > 0.) 

One notices that when n = 0, the limit is independent of the expected background b. 
This is required in the Bayesian approach: we know that exactly zero background events 
were produced (when no events at all were produced), and this knowledge of what did 
happen makes what might have happened superfluous. An interesting corollary is, in 
the case of no events observed, uncertainties in estimating the background rate are 
of no consequence in the Bayesian approach, and must not contribute any systematic 
uncertainty to the limit. This reasoning does not hold in the frequentist framework, 
where what might have happened definitely does influence the limit. 

For comparison, limits for fixed e = 1 with 6 = or 6 = 3 are also shown in Tabled 
It is interesting that these two columns start out equal at n = and differ by almost 
exactly 3 for n > 11. In contrast, the difference between the b = and 6 = 3 columns 
for e = 1.0 ±0.1 is already greater than 3 at n = 6, and continues to grow as n 
increases; it is not clear whether the difference approaches a finite value as ^ oo. In 
any case, the limits for e = 1 exactly are all smaller than the corresponding limits for 
e = 1.0 ± 0.1, as expected. 

4.3 Coverage 

The main quantity of interest in this subsection is the frequentist coverage probability 
C as a function of Strue (for fixed etrue and b) . Because both the main and the subsidiary 
measurements involve observing a discrete number of events, the function C(strue) will 
have many discontinuities. On the other hand, C(etruc) will be continuous (for fixed 
Strue)- The explanation of this effect is as follows: 

The measured data are n events in the main measurement and m events in the 
subsidiary measurement. For each observed outcome (n,m) there is a hmit s^{n,m). 
This hmit includes the effect of marginahzation over e. 

An (n, m) with n > and m > are possible, and the probability P of observing 
(n, m) can be calculated as the product of two Poissons. (It will depend on struo, Ctrue, 
. . . ) If we look at all the possible limits we can obtain, 

{s^{n,m)\n > and m > 0} (19) 

and sort them in increasing Su, the are countably infinite in number and dense in the 
same way that rational numbers are dense in the reals. 

To compute the coverage as a function of Struej we simply add up all the probabilities 
of obtaining (n,m) with Su(n, m) > Strue: 

C ~ P{n,m) A — {(n, m)|strue < Su{n,m) and 7i > and m > 0} (20) 

This sum is over a countably infinite number of terms. If we increase Stmc slightly to 
Struc + ds and recalculate the coverage, we have to drop all the terms 

{{n, m)|struc < Su(n, m) < Struo + ds} (21) 
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from the previous sum (the P(n,m) for each term also changes continuously with Struej 
but this is no problem). If there are M > such terms, there are M discontinuities in 
the coverage in the interval [strucStrue + ds], since P{n,m) for each of these is finite, 
and we lose them one by one as we sweep across the interval [struei Struc + ds]. 

But it seems that, in general, we can always find a solution to Struo < Su{n,m) < 
Struc + ds for finite ds by going out to larger and larger n and m. So, although the 
discontinuity may be tiny, we can always find a finite discontinuity in any finite interval 



On the other hand, if we keep Struo fixed and vary etruo, we always sum over the same 
set of (n,m), since the definition of A does not involve Ctruc, and P{n,m) is continuous 
in etruc- So the coverage is continuous as a function of etruc for struc fixed. 

Plotting a curve that is discontinuous at every point is somewhat problematicaL 
The solution adopted here is to plot the coverage as straight line segments between 
the discontinuities, ignoring any discontinuities with |AC| < 10~^. Figure El shows 
C(strue) for the case = 90%, etrue = 1, nominal 10% uncertainty of the subsidiary 
measurement of e, and 6 = 3. We observe that C(strue) > (3 m. this range, and it is 
not clear numerically whether C(strue) — > /5 as Struc ^ oo. The same conclusions hold 
for Fig. ini which illustrates the same situation with a 20% nominal uncertainty for the 
e-measurement . 

Figure m shows C(etrue) for /3 = 90%, Stme = 10, k = 100, and h = 3 — continuous 
as advertised. The shape of the curve is quite similar to that of Figs. El and El so it 
seems that the coverage probability (with b fixed) is approximately a function of just 
the product of etrue and Strue- This approximate rule is likely to fail in the limit as 
Ctruc and Struc oo, for example, but it seems to hold when etrue and Strue are at 
least of the same order of magnitude. 

When etrueStrue is small, of order 1 or less, the coverage is ~100%, as in the simple 
case of Fig. [H Otherwise, the behavior of coverage in Figs. EHZl is superior to that of 
Fig. H which has a much larger amplitude of oscillation. 

Another frequentist quantity that characterizes the performance of a limit scheme 
is the sensitivity, defined as the mean of s^. Figure |H1 shows the sensitivity as a function 
of Strue for the case of Fig. El (su) is observed to be nearly linearly dependent on Strue- 
There is one complication here: when the subsidiary measurement observes m = 
events, and the prior for s is fiat, = oc. Since the Poisson probability of obtaining 
m = is always finite, (su) is consequently infinite. So we must exclude the m = 
case from the definition of (su). (In Fig. El the probability of obtaining m = is 



4.4 Other priors for s 

A weakness of the Bayesian approach is that there is no universally accepted method 
to obtain a unique "non-informative" or "objective" prior p.d.f. Reference jO], for 
example, states: 

Put bluntly: data cannot ever speak entirely for themselves; every prior 
specification has some informative posterior or predictive implications; and 
"vague" is itself much too vague an idea to be useful. There is no "objective" 
prior that represents ignorance. 



of St, 



;ruc- 
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Nevertheless, Ref. jH] does derive a 1/a/s "reference prior" for the simple Poisson case, 
which is claimed to have "a minimal effect, relative to the data, on the final inference". 
This is to be considered a default option when there are insufficient resources for 
detailed elicitation of actual prior knowledge" . 

Reference [10^ attempts to discover the optimal form for prior ignorance by con- 
sidering the behavior of the prior under reparameterizations. For the case in question, 
the form 1/s clearly has the best properties in this respect. 

We are using an flat (s") prior for this study, which seems to be the most popular 
choice, but the Appendix works out the form of the posterior using an prior, so 
we can briefly here summarize the results for the 1/s and cases: 

The 1/s prior leads to an unnormalizable posterior for all observed n when 6 > 0. 
The posterior becomes a 5-function at s = 0, Su = for any j3, and the coverage is 
consequently zero for all Strue > 0. This clearly is a disaster. 

The I/a/s prior results in a posterior p.d.f. qualitatively similar in shape to those 
of Fig|l[ except that the p.d.f. is always infinite at s = 0. For n ^ 6, this produces an 
extremely thin "spike" at s = 0, which has a negligible contribution to the integral of 
the posterior p.d.f. A more significant difference (for frequentists) between the l/^/s 
and the s° case is that the coverage probability is significantly reduced: for the case 
of Fig. El the 1/a/s prior pushes the minimum coverage down to ~0.87. So the l/y/s 
prior leads to violation of the frequentist coverage requirement; it undercovers for some 

values of Strue- 

One might also seek to further improve the coverage properties by adopting an 
intermediate prior. For example, an prior would reduce the level of overcoverage 

obtained with the s° prior. How acceptable this approach would be within the Bayesian 
Statistical community is an interesting question. 

It should be noted that all the prior p.d.f. 's considered in this note are "improper 
priors" — they cannot be correctly normalized: In the case of the s° and l/i/i priors, 
the integral from to any value sq is finite, while the integral from sq to infinity is 
infinite. The corresponding integrals of the 1/s prior are infinite on both sides for all 
So > 0. Improper priors are dangerous but often useful; "improper posteriors" are 
generally pathological. Extra care must be taken when employing improper priors to 
verify the normalizability of the resulting posterior — when using a numerical method 
to obtain the posterior, it is very easy to miss the fact that its integral is infinite. 

4.5 Restrictions 

We summarize here the restrictions forced on the priors for s and e — see the Appendix 
for the analytical causes. The discussion below assumes 6 > 0. The prior for s being 
of the form s"~^, we must require a > 0, as discussed above. 

As specified in this note, the prior for e, being taken from the posterior from the 
subsidiary measurement with a fiat prior, has been given no freedom. Should the 
subsidiary measurement observe m = events, the posterior for s is not normalizable 
when a > 1: Su = c>o when m = and a > 1. 

This behavior is due to a well known effect: the e prior becomes Ke~'^'^ when m = 0, 
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which remains finite as e ^ 0. All such cases^^ yield s^^ = oo when a > 1; any positive 
a < 1 cuts off the posterior at large s sufficiently rapidly to render it normalizable. 
From this point of view, a l/^/s prior may seem preferable, but on the other hand, 
having Su = oo when m = seems intuitively reasonable. (In general, we have s,^ = oo 
for m < a — 1, but a >2 are not popular choices.) 

There is another approach possible to the gamma prior for e: we may simply specify 
by fiat the form of the prior as 

Pi^\f^)de = -^-Yij;^ de (22) 

where fi is no longer required to be an integer. In practice, one then might obtain 
and K from a subsidiary measurement whose result is approximated by the gamma 
distribution. In such cases, one must require /i > a to keep the posterior normalizable. 
Note that in this form, fi/n is the mean of the e prior, [ij, — 1)/k is the mode, and fi/K^ 
is the variance. The subsidiary measurement is often analysed by other experimenters, 
who chose statistics to quote for their central value and uncertainty (omitting additional 
likelihood information). It is then important to obtain /i and k in a consistent way 
from the information supplied by the subsidiary measurement. If e, for example, were 
estimated by a maximum likelihood method, one would identify the estimate with 
(/i — 1)/k rather than fi/n. 



5 Conclusions 

Results have been presented on the performance of a purely Bayesian approach to the 
issue of setting upper limits on the rate of a process, when n events have been observed 
in a situation where the expected background is b and where the efficiency/ acceptance 
factor e ± (Te has been determined in a subsidiary experiment. We find that this 
approach, when using a fiat prior for the rate, results in modest overcoverage. Plots of 
the expected sensitivity of such a measurement and of the coverage of the upper limits 
are given. It will be interesting to compare these with the corresponding plots for other 
methods of extracting upper limits, to be given in future notes. Reference provides 
the limit calculating software associated with this study in the form of C functions. 



6 Appendix A — Analytical Details 

Here we present the details of the analytical calculation of the posterior p.d.f. for s. 
For generality, we work through the calculation with a prior; a fiat prior is then 
the special case a = 1. 



6.1 Posterior for s with e and b fixed 

We measure n events from a process with Poisson rate es + b, and we want the Bayesian 
posterior for s, given improper prior We compute the posterior for fixed e and b 

^^A Gaussian truncated at e = is the standard example. 
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in this subsection; the calculation with our prior for e follows in the next subsection. 
We have ^ 

posterior: p(s|e, 6, n)ds = — re~"(es + b)^s°'~^ds 

where all factors not depending on s have already been absorbed into the normalization 
constant Ms, which is defined by 

poo Un+a poo 

Us= e-'\es + hYs''-^ds= e-^''u''-^{l + uYdu {u = se/b) 

Jo Jo 

where we have performed the indicated change of variable. 

Expanding (1 + m)"" in powers of u using the binomial theorem, we get 

^ ^ "^(n-ky.kl ^ (n-k)\k\ 

k=0 ^ ' k=0 ^ ' 



Recognizing this as of the general hypergeometric form, we write it as 
Ms = e-"r(a + n) 



b nin — 1) b"^ 



n ^ , 

^ a + n - 1 1! ^ (a + n-l)(a + n-2)2! ^ ' 



to make the hypergeometric nature more explicit. Using the modern notation^21 for 
the falling factorial 



this is expressed as 

n k i,k 

Ms = e-Tia + y — --— = e-"r(a + n)M(-n, l-n-a,b) 

^-^ (a + n — 1)- /c! 

k=o ^ ' 

where M is the notation of • iM-, a confluent hypergeometric function, is often 
written \F\, and the relation given here is only valid for integer n > 0.) Note that 
M{—n, 1 — n — a,b) is a polynomial of order n in 6 (for n a non-negative integer), and 
is related to the Laguerre polynomials. When a = 1, we get M{—n, —n,b), which is 
related to the Incomplete Gamma Function. When a = 0, we get M{—n, 1 — n,b), 
which is infinite, so we require that a > 0. 

Our posterior probability density for fixed e is then given by 

p{s e, 6, n)ds = — — ■ — — — -ds 

i (a + njM[—n, 1 — n — a,bj 

6.2 Posterior for e of the subsidiary measurement 

The subsidiary measurement observes an integer number of events m, Poisson dis- 
tributed as: 

P(m\e) 
^ ' ^ ml 
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where k is a real number (connecting the subsidiary measurement to the main measure- 
ment) whose uncertainty is neghgible, so k can safely be treated as a fixed constant, k 
might be thought of, for example, as based on a cross section that is exactly calculable 
by theory. There is negligible (i.e. zero) background in the subsidiary measurement. 
The prior for e is specified to be fiat. The Bayesian posterior p.d.f. for e is then 

p{e\m) 



(or r(m + 1) instead of m! in the denominator if you prefer). This is known as a gamma 
distribution. 

The mean and rms of this posterior p.d.f. summarize the result of the subsidiary 
measurement as: 

m + 1 , \/m + 1 

e = ± = eo ± (Te 

K K 

Note that the observed data quantity in the subsidiary measurement is an integer m, 
while the quantity being measured by the subsidiary measurement is a positive real 
number e. 

6.3 Posterior for s with gamma prior for e (6 fixed) 

Next we compute the joint posterior p(s, e|6, n)dsde using the s"~^ prior for s and our 
gamma distribution prior (i.e. the posterior derived from the subsidiary measurement) 
for e > 

r / \ 7 [Kdj^C dc -1 / / \ 2 / 2 

prior for e: 7r[e)de — — — /j, — m + 1 — (eo/c"ej K' — eo/C( 

i [fj,) e 



€ 



where it is convenient to write for m + 1. We have for the joint posterior p.d.f. 
p{s, e\b, n)dsde — — ^7r(e)e~^*(es + bYs"'~^dsde 

•1^ s,e 

where 

/»oo /*oo /»oo 

Ns,e= / / 7r(e)e-^^(es + 6)"s"-^ciscie= / n{e)J\fsde 
Jo Jo Jo 

We calculated jV^ above, so we have 

poo 

Ns,e = r(Q; + n)M{-n, l-n-a,b) / e~"7r(e)de 

Jo 

Ms,e = «'T(// - a)T{a + n)M{-n, 1 - n - a, 6) /r(^) 

WIS, e 0. n)dsde = — — — — — —dsde 

' ' ^ r{ii-a)T{a + n)M{-n,l-n-a,b) 

The marginalized posterior for s can then be expressed as 



p{s\b, n)ds 



oo 



p(s, e|6, n)de 







r(// - a)V{a + n)M{-n, l-n-a,b) 
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where the integral is given by 

Jo 

The same procedure that was used for the normahzation integral can be applied here, 
producing 

Jo 



s"n! ^r(/i + n-A:) 



hn 



{s + Ky+'^^ {n-k)\k\ 



b{s + k) 



k 



^ (g + ^)/.+n ^(/^ + ^Wi-n, l-n-i^,bis + k)/s) 

, T(a + n) s'^+^'-^k''-" M(-n,l-n- a.b(s + k)/s) , 

p{s\b,n)ds = —, r^r- — — ; '■ ds 

' ' ^ T{fi- a)T{a + n) {s + Ky+'' M{-n,l - n - a,b) 

which has a particularly simple form when the background term is zero: 

pis = 0, njds — — — — r — ; — ds 

' ' ^ r(/i-a)r(a + n) (s + 

a Beta distribution of the 2nd kind. Note that we must require /i > a > to obtain a 
normalizable posterior. 

Our posterior p.d.f. for s with e (and b) fixed is recovered exactly by taking the 
limit of p{s\b,n) as (7g 0. This means that the limit of as —>■ is identical to 
the value of Su when e is known exactly. This property may seem obvious, but it is 
violated by some frequentist methods of setting limits, so it is worth mentioning. 

6.4 Calculating the limit 

Wc need to integrate p{s\b,n) up to some limit Su, which can be done analytically as 
follows. 

Jo V{ii-a)T{a + n) Jo ^ ' M{-n,l - n - a,b) 

where the substitution t = has been performed. Re-expanding the polynomial M 
and integrating term by term yields 



1^ Ix{a + n — k, II — a)n-b'' / v-^ n- b^ 

— / k — 

where is the standard notation for the Incomplete Beta Function 
which also satisfies the following recursion: 



X 



16 



6.5 Integer moments of the marginalized posterior 

Using the same technique as above, we can calculate the jth moment of the posterior 
p.d.f. as 

{s-^} = / s-^p[s\b, n)ds 



{fi-a-iy- M{-n,l-n-a,b) 
where we utilize the rising factorial notation [T^ 

The expression for the mean of the posterior when a = 1 can be simplified using 
the identity 

/ b \ 6"+^ 
M(-n, -n - 1, 6) = 1 M(-n, -n, b) + 



n + lj ' ' ' ' {n + l)\ 



obtaining 



/t(n+l-6) , Kb-+' 
mean(a = 1) = (s)|a=i = h 



/i-2 (/i-2)n!M(-n, -n,6) 



Note that the 2nd term is very small when n ^ b. 
The recurrence relation lil? 



r(r — l)M{q, r — l,z) + r(l — r — z)M{q, r, z) + z{r — q)M{q, r + l,z) = 
leads to a recurrence relation between moments 

{s^) = -i^ + n+j-l-b) ^ n%{a+j-2) 

12- a- j {/J,- a- j + l){fj.- a- j) 

The special case a = 1 then yields 

'2 r (2 + n - 6)6"+^ 



(/i-2)(/i-3) 



{2 + n-b){l + n-b) + b + 



n\M{—n, —n, b) 



which leads to this approximation for the variance of the posterior 

K^il + n) K^{l+n-bf . 

variance(a = 1) ^ (^_2)(^-3) + - 2^/. - 3) 

6.6 Posterior for s with gamma priors for e and 6 

Here we very briefly consider the case where the background parameter b also acquires 
an uncertainty. This case is more general than the fixed b case that is the main subject 
of this note: the fixed b case will be the subject of additional studies employing various 
popular frequentist techniques, with the goal of comparing their performance. We 
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judge the more general case considered in this subsection to be more comphcated than 
necessary for the purpose of comparing the various methods, but it is instructive to 
document the fact that the Bayesian method can easily handle the more general case. 

We assume a 2nd subsidiary measurement observing r events (Poisson, as was the 
case for e), which, when combined with a flat prior for b, results in a gamma posterior 
for b of the form 

t\ 

where is a calibration constant (analogous to k in the subsidiary measurement for e) . 

The posterior for b becomes the prior for b in the measurement of s. After deter- 
mining the joint posterior p(s, e, b\n) by using our priors for s, e and b, we marginalize 
with respect to e and b, resulting in 

,,,, T(a + n) s"+"-iK^-"F(-n,p;l-n-u;(s + K)/(sc<;)) , 

P[S\n)dS = — r ; — ; -— dS 

T{n — a)T{a + n) (s + F{—n,p;l—n — a;l/uj) 

where we write p = r + 1 for convenience, and F is the hypergeometric function |14j. 
As long as n is a non-negative integer and a > 0, F{—n, p; 1 — n — a; x) is a polynomial 
of order n in a; (closely related to Jacobi polynomials). 

This marginalized posterior for s can then be integrated, with the result 

Z"**" Ix{a + n - k, jj, - a)n-p'' uj~^ / ""■"/'^ f \ 

Jo ^^''""^ ^ («+n-i)^ / ^ (a+n-i)^Tr = t:^) 

These two equations closely resemble the main results of sections 16.31 and 16.41 to 
recover the fixed b results, simply substitute buj for p above, and take the limit u oo. 



7 Appendix B — Average Coverage Theorem 

In this appendix we prove that Bayesian credible intervals have average frequentist 
coverage, where the average is calculated with respect to the prior density. We start 
from the Bayesian posterior density: 

/ I X Pin I s) tt(s) , , 

V{s\n) = ) \ W . (23) 

Jg F[n I s) 7r[s) ds 

For a given observed value of n, a credibility-/? Bayesian interval for s is any interval 
[shin), s-[j{n)] that encloses a fraction P of the total area under the posterior density. 
Such an interval must therefore satisfy: 

(3 = / p{s\n)ds, (24) 



or, using the definition of the posterior density: 



su{n) POO 

P{n\s)7r{s)ds = P Pin\ s) 7r{s) ds. (25) 

L(n) Jo 
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Now for coverage. Given a true value St of s, the coverage C{st) of s\j{n)] is the 

frequentist probabihty that St is included in that interval. We can write this as: 

C{st) = Yl Pi^\st). (26) 

all n such that: 
SL{n)<s<su{n) 

Next we calculate the average coverage C, weighted by the prior 7r(s): 

POO 

C = / Cis)7r{s)ds, 







P(n I s) 7r(s) (is, using equation (j26|l . 

- all iv) oii/^Vi +Vici + - 



all n such that 
SL(n.)<s<su(") 



°° rs-[j(n) 

/ P(n I s) 7r(s) (is, interchanging integral and sum 

n=0 >^«L(ri) 
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/"OO 

n=0 



P{n I s) 7r(s) (is, using equation 



/'OO ^ 

/3 / P(ri I s) 7r (s) (is, interchanging sum and integral, 

n=0 



oo 



/3 / 7r(s) (is, by the normalization of P(n I s). 





= /3, by the normalization of vr(s). 

This completes the proof. We have assumed here that the prior tt{s) is proper and 
normalized to 1, but the proof can be generalized to improper priors such as those we 
considered in this note. A constant prior for example, can be regarded as the limit for 
Smax — > C)0 of the propcr prior: 

Vr(s I Smax) = , (27) 

"5 max 

where '(9(x) is if x < and 1 otherwise. We then define the average coverage for the 
constant prior as the limit: 

POO 

C = lim / C(s) 7r(s | s^^ax) ds. (28) 

^max ^ -t-OO I 



"'^^The best way to understand this step is to draw a diagram of s versus n: one is integrating and 
summing over the area between the curves sl(".) and su{n). The hmits on the sum and integral 
depend on the order in which one docs these operations and can be derived from the diagram. 
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The previous proof can now be applied to the argument of the hmit and leads to the 
same result. 

The average coverage theorem remains valid when s is multidimensional, for exam- 
ple when it consists of a parameter of interest and one or more nuisance parameters. 
In that case one needs to average the coverage over all the parameters. 
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Stmt 



Figure 1: Coverage as a function of the true signal rate s for Bayes 90% limits, for the 
simple case of no background and no uncertainty on e = 1. The dotted line at C = 0.9 
is given to show that the coverage never falls below 90% (in this simple case). 
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Figure 2: Comparison of our discrete probability for eo (shown as a histogram, see 
eqn. (fTH|l ) and Gaussian (continuous curve) for the case e = 1.0 ± 0.1. 
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Figure 3: Comparison of our likelihood (dashed, see eqn. (|TT)|) ) and Gaussian (solid) 
for the case e = 1.0 ± 0.1. 
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Figure 4: Posterior densities p{s\b,n) vs s for n — 1, 3, 5, 10. In each case, b — 3 and 
e = 1.0 ± 0.1 (i.e. k = 100 and m=99). 
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Figure 5: Coverage of 90% upper limits as a function of Strue for etrue = 1) nominal 10% 
uncertainty of the subsidiary measurement of e, and expected background 6 = 3. 
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Figure 6: Coverage of 90% upper limits as a function of Strue for etrue = 1) nominal 20% 
uncertainty of the subsidiary measurement of e, and expected background 6 = 3. 
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Figure 7: Coverage of 90% upper limits as a function of etrue for Strue 10, nominal 
10% uncertainty of the subsidiary measurement of e, and expected background 6 = 3. 
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Figure 8: Sensitivity of 90% upper limits as a function of Strue for ^true — 1; nominal 
10% uncertainty of the subsidiary measurement of e, and expected background 6 = 3. 
For reference, the sensitivity for Ug = is also given (dashed). 
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